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1. Introduction 

The gg — ► Z°/^*bb — > //£>& process is of high experimental interest in view of the forthcom- 
ing LHC experiments, since it e.g. represents an irreducible background to the 'gold-plated' 
Higgs channel H — > ZZ^ — ► 4 leptons (see e.g. [0). Historically, the process has been 
fully calculated at tree level by Kleiss and Stirling Q and for the first time successfully 
implemented inside the AcerMC || Monte-Carlo generator. In the present time there is 
a plethora of other Monte-Carlo generators implementing this process in various advanced 
fashions (see e.g. [Q]), however a few issues still need to be resolved in a consistent fashion: 

• At tree level the gg — ► Z°/j*bb — ► ffbb process is actually NNLO in a s with respect 
to the order a® 'pure' Drell-Yan process bb — > Z°/j* — > // . The same final state 
can thus be achieved by using the Sudakov parton showering procedure, which by 
definition re-sums the large logs of the order a s In M| / which burden the higher 
order corrections as the gg — ► Z°/^/*bb and might thus be better at least in the 
'intermediate' kinematic regions (see e.g. || where the gg — ► Z°/7*66 — ► //66 process 
was actually removed with this argument). While the argument by all means stands 
it would nevertheless be preferable to have the processes consistently combined in 
orders a", n = 0,1, 2, achieving the undisputed validity over the whole phase-space. 

• The b-quarks are reasonably massive; still it is customary to treat all partons in- 
coming to the hard process as massless, which is strictly speaking consistent only if 



- 1 - 



also the final state b-quarks are also treated as massless. In 'full' Monte-Carlo proce- 
dures, where the produced b-quarks are further hadronized into jets, the mass of the 
b-quark needs to be present and is thus added in and ad-hoc fashion. Furthermore, 
neglecting the heavy quark masses of the incoming quarks has been shown to have 
an observable impact and can in fact be consistently added into the Factorization 
Theorem g |, |. 

The existing prescriptions deal either with massive particles |Q, g", P on the level of in- 
tegrated cross-sections or with explicit Monte-Carlo algorithms involving light (~ massless) 
particles (e.g. ||, [l(], 11 , [l^, [D| ) while a first attempt at the combination of the two was 



attempted in the paper [14], where an algorithm combining the two features was developed 



but implemented only in terms of order a] correction while for the gg — > Z°/'y*bb process 
the order a 2 combination procedure is needed. The aim of this paper is to show that the 
procedure developed in [14] is in fact iterative and can thus provide a consistent proce- 
dure applicable (at least) at tree-level which can be and is implemented in terms of a full 
gg — > Z /j*bb Monte-Carlo procedure. 



2. Combining the Perturbative QCD and Sudakov Showering in Massive 
Hadroproduction 

2.1 Theoretical Basis 

Let us start by repeating the considerations of the 
Factorization Theorem [||, 15, [h], 17, 18 1, which states that 
the hadronic cross cross section ctab^x is related to the 
perturbatively calculated (e.g. by using Feynman diagram 
technique) parton level-cross section cj a b^x by: 



a,b 



fa/A ® &ab->X ® fb/B, 



(2.1) 




f./B 



with /j/j = f i j I {x,^F) denoting the parton density func- 
tions (PDF), giving a probability that a fully evolved par- 
ton i is produced by the parent hadron I, at the factoriza- 
tion scale [If with a certain energy fraction x. The PDFs 
are convoluted with the hard parton-level cross section 
a a b^,x which can in general differ from the perturbative 
(pQCD) parton level-cross section cr a b^x- 

Note that the term hard cross section in the above application of the Factorization 



Figure 1: Schematic represen- 
tation of the example of a gluon 
splitting g — > H H combined with 
a hard subprocess a . 



Theorem demands that the hard cross section expression of Eq. 2.1 is indeed describing 
the hard ('short-distance', high energy) process only, i.e. all the 'long-distance' contribu- 
tions in form of collinear /mass singularities and corresponding large logarithms in form of 
a s log(fi F /m 2 ) need have been explicitly subtracted since they are already included (re- 
summed) in the PDFs |], 0, g]. The factorization scale /ip sets the dividing limit between 
the two kinematic regimes. 
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Let us illustrate how this applies to the case studied in this paper: In a perturbative 
calculation of order a" an incoming gluon splits to a heavy quark pair g — > HH; the other 
incoming parton we mark as the spectator s (c.f. Figure ||). The Factorization theorem 
then gives: 



(7 



• AB^XH - hi A ® ^ s i xS ® /s/B ■ (2-2) 

Possible summations and permutations of incoming flavors are omitted and only the con- 
volution with the parton density functions remains explicit. 

Stipulating that the (hard/soft) scale /j, is set by the heavy quark kinematics (e.g. heavy 
quark propagator; other choices are possible) then : 



If the scale is hard enough, fi > \ip , the pQCD calculation need not be modified, i.e. 

(«) (n) 
gs^XH gs^XH 

in this kinematic region. 



si!** = ^ixu o* > mf) (2.3) 



• if the scale is soft, {j, < [ip, the heavy quark production in gluon splitting is 'long- 
distance' and thus already included in the PDF fjj/i- In other words, the incoming 
gluon cannot be resolved at this scale [ip and the heavy quark H should be considered 
as the incoming (fully evolved) parton. One should thus use an alternative calculation 
^Hs-^x w ith an incoming quark H instead: 

o^AB-ix = fn/A ® ^Hs^X ® fa/B (A* < MpO (2-4) 
and correct the perturbative calculation correspondingly. 

The latter case (fi < fip) however needs to be explained a bit further: One should 
not simply set o^l^xg = since only the collinear kinematic limit (explicit large logs 
a s log(/i^/m 2 )) is resummed in the fn/i, whereas such logs are only a part (limit) of the 

* (n) 

full pQCD calculation & s _^ x jj and it is only these contributions that need to be removed 

in order to prevent double counting of the inclusive Ojj s 2}x contribution in this kinematic 
region. Instead one should put: 

h (n) _ On) _ On) , 2 

which can then be used in the hadronic cross-section expression: 

CT AlLxH = hi A ® ^gixfl ® = /ffM ® ^gixfl ® " hi A ® ^ubt ® /»/B • ( 2 -6) 



One can immediately deduce one property of the subtraction term from Equation 2.3: 



^bt = (M>MF). (2.7) 

The full (correct) cross section involving the heavy flavour excitation in the initial 
state (category HE in the paper ||) thus needs the addition of the a^ s ^ x hard process to 
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Figure 2: Schematic representation of heavy quark H entering the perturbative calculation at 
order (n-1) as fully evolved (left) or participating internally via the splitting of the incoming gluon 
at order (n) (middle) with the appropriate subtraction term (right). The subscript s denotes the 
other incoming parton (which can also be a heavy quark) and X the inclusive final state. 



correctly cover the soft region: 

CT AB->X(H) = fg/A ® *£l»xH ® f»/B ( 2 - 8 ) 
" fg/A ® °^lt ® fs/B 

+ fu/A ® O-Hs^X ® fs/B , (P < PF) 

as illustrated also in Figure ||. 

At this point it should be emphasised that all the above hadronic cross-sections are 
effectively of the same order if one considers that the heavy quark density f H /j is of the 
effective order a s higher with respect to the dominant (gluon or valence quark) parton 
distribution functions (see e.g. || for details). Also note that the final state heavy quark 
H need not be resolved in the final state since both final states XH and X from the two 



hard contributions participate in the full expression of Equation |2.8| and the H can in the 
second case appear only in the soft (parton showering) processes. 

An excellent basis for deriving the desired rules and expressions is the paper of Olness, 
Scalise and Tung ||, which presents a comprehensive review of the formalism needed for 
obtaining the fully infra-red safe hard cross-sections with the emphasis on isolating and 
subtracting the divergences related to the heavy quark terms. The derivation of || (or 
equivalently and our previous paper [|l4|]) gives the expression for the subtraction term 
of the order (n) as (please consult the Appendix ^ for details): 

= (/*</*) (2.9) 

where the first-order in a s term /jyL is the perturbatively calculated parton distribution 
function fty of a parton inside another parton which is explicitly given as (c.f. |l4|) : 

<«<>) 

with Pj^i being the well-known splitting kernels and also contains the explicit massive 
divergence logarithm. 
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This is not necessarily the final step one needs to make to be able to evaluate the 
hadronic cross sections, since e.g. the spectator s can also be a gluon splitting to heavy 
quarks (as is indeed the case in the studied gg -> Z°/^*bb process) and thus the whole 
procedure and calculation of the subtraction terms needs to be recursively (iteratively) 
repeated for both n th and (n-l) th order terms in order to obtain all the missing subtraction 
terms. The repeated procedure is somewhat lengthy and is thus presented in the Appendix 
|A| of this paper. 

Note that the subtraction term from Equations |2.8| , |2.9| can be combined with either of 
the two hard cross-sections; joining it with (n-l) th (Eq. |J) order cross-section instead of 
n th (Eq. ^6|) one gets: 

^AB— >X = Uh/A ~ fg/A ® fgj g ) ® ^Hs^X ® fs/B = Jh/A ® ^Hs^X ® fs/B , i. 2 - 11 ) 

where the new parton distribution function Jh/A = {fH/A~fg/A®f ( H > / g ) nicely expresses the 
physical origin of the subtraction term. In addition to being the collinear limit of the n th 
order calculation, it also represents the first (fixed) order component of the QCD evolved 
parton distribution functions which is thus explicitly removed from the fully re-summed 
(to all orders in a s ln(/j, 2 /m 2 H )) function fu/A- 

The experimentalists are generally interested in the fully differential cross-sections from 
which exclusive topologies (events) can be picked. The procedure to obtain an exclusive 
topology from an inclusive differential cross-section is well established in various flavors 
HI, H|, and is commonly known as the Sudakov parton showering. Using the case of Eq. 
|2.4j as an example, the Factorization Theorem states that the incoming heavy quark is fully 
resolved at the factorization scale \i = (jlf- The probability dP g ^ H jj(n) of the heavy quark 
un-resolving back to the gluon (back-evolution) at a lower scale (j, < /J,p (and producing 
an exclusive state with an additional H) is then given by the differential of the Sudakov 
exponent (see e.g. ||): 

For the given example of gluon splitting to heavy quarks with dP g ^ H g(fj>) = dSg-tHifJ*) 
one gets: 

d <™ Ib = dS g -+ H (v) Sh/aM da^l /. /b (mf) , (2-13) 

which is also written as fully differential in all variables. Note that the dS g ^u{n) term 
in the Monte-Carlo context represents an explicit showering step (i.e. a parton-showered 
parton level 'hard' event). A fact also worth noting is that this expression is now of the 
order n and contains one resolved particle H, thus giving the same configuration as the 
perturbative n th order expression of Eq. |2.6| in the region (// < /xp). There is clearly an 
overlap between the results of the two in this region, producing the 'double-counting' in the 
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overlap region. The subtraction term of Eq. however retains its role also in differential 
form: 

d^l = (M < ) (2.14) 

and this fully differential subtraction term depends on the variable \x through: 

dz d/u, 2 



(2.15) 



where the variable needs to be explicitly kinematically related to the mass singularity 
in the (differential) n th order cross-section. The Equation 2J3 can now be rewritten in 
differential form: 



da 



(n) 

AB^XH 



fg/A ^iixH f*/B ~ fg/A daK 'sLxB f S /B - fg/A da sult fs/B ■ 



(n) 



» 



'gs-XH 



gs^XH ■ 



(2.16) 



Correspondingly, anticipating the exclusive final state containing one heavy quark, the 
differential form of Equation |2.8| is then given by: 



da 



AB^XH 



fg/A dff^xfl 



fs/1 



(2.17) 



- fg/A da^l t f s/B 

I j shower 

AB^XH ' 



(A* < V-f) 



which, when inserting the explicit expressions from Equations 2.13 and 2.14 gives: 

' in \ m fs/BM (2.18) 

<•(!) 



d^AB^XU = fg/A(^F) da^ xn f s/B (fJ, F ) 



fg/AiVF) tin/git*) ^HsJx fs/BifJ-F) 



(/i < //f) 



+ dS g ^ H (v) Ih/a^f) ^hs^x /s/b(mf) • (m < Mf) 

It is trivially obvious that the subtraction term still cancels the collinear limit (mass 
divergence) of the perturbative n th order calculation of Eq. 2.16| , while its impact on the 
(n-l) th order (showered) calculation is again best seen by combining the last two lines of 
Equation 2.18| and results of Eqns. 2. IS and 2.14 into: 



da 



(n-l) 
AB^XH 



dSg^ H (^) fH/A(^F) ~ fg/A^F) <VH/g(p) 



da 



(n-l) 
Hs^X 



/s/s(w 



Taking the limit /i — > fip one quickly sees that: 



dS g -+n{}l) fH/AiPF) M > f^F fg/A{^F) aS ^ F Pg-+H d ® 
fg/AifiF) df§J g in) fi-> (IF fg/AiPF) "'^ Pg->H d&, 



(2.19) 

(2.20) 
(2.21) 



with d& denoting all the variables in the differential. In this limit [i — > jxp the two 
terms thus cancel on paper, i.e. exactly. Since \ip is by definition the highest virtuality 
reachable by parton showering approach, the kinematic region fi > fj,p is thus populated 
solely by the n th order contribution and the continuation on the transition point is smooth. 
Choosing a kinematic relation which defines fi in terms of kinematic quantities measuring 
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the collinearity (e.g. transverse momentum of the splitting, virtuality of the participating 
particle/propagator) one achieves an ordered subtraction scheme ranging from the collinear 
region to the showering limit \xp. 

The result of these deliberations is that the subtraction term besides it's obvious role 
actually smoothly interpolates between the parton-shower evolved (n-l) th order and n th 
order perturbative expressions while removing the double-counting contributions in the 
overlap region (/i < //p). 

There is no obvious reason why this procedure cannot be made iterative since the 
subtraction terms are designed to cancel all mass divergences present in a perturbative 
calculation of a certain order. In fact, a consistent iterative (and recursive) procedure has 



been developed and implemented in our paper [14], however the iterative nature of it was 
not tested in the examples developed thus far. The process of the associated Z boson 
production with two heavy (b) quarks gg — ► Z° / / ybb serves as a good test case when one 
wants both b-quarks to be resolved (observable) and correctly described over the whole 
phase space. 

2.2 The Implementation of the Showering and Overlap Removal 

In order to implement the above procedure in a Monte Carlo algorithm one needs to define 
a specific mapping of the variables fj,, z in the Sudakov (parton) showering algorithm (c.f. 



Equation 2.13 ) as well as appropriate kinematic transforms between the four- momenta of 
the partons undergoing the showering. 

The choice of the kinematic mappings and transforms to suit our procedure is in 
principle arbitrary, as long as it properly takes into account the heavy parton masses. 
Consequently the defined procedure could be matched to the parton showers of e.g. Pythia 
|H or Herwig [19] if the mass of the partons were incorporated. 

Nevertheless, we chose to implement our own showering algorithm based on the one 
suggested by Collins et al [11, p^] and extended it to properly incorporate heavy 



parton masses. The explicit derivation of the showering algorithm, kinematic mappings 
and transforms applied to our procedure have already been described in our previous paper 



1 14]. Summarized briefly, the implemented showering algorithm has the evolution variable 
H chosen to be equal to the virtuality of the particle (in our case heavy quark) and in each 
showering step the rapidity of the showered system is preserved, which sets the variable z 



of Equation 2.12. The value of the factorization scale //p is not pre-determined (i.e. can be 
picked from one of reasonable options, e.g. the invariant mass of the subsystem or similar). 

The choice of the showering algorithm implementation was motivated by two main 
points: 



This subtraction term of Equation 2.14, albeit derived in a different way, is in form 
identical to the subtraction term derived by Collins et al. This fact motivated us 
to implement the parton showering algorithm according to prescriptions given in the 
cited papers. 

The selected kinematic setup is motivated by the fact that, as shown by Collins et al 
[ |Tl| , |20| , pi , 22], the procedure of parton showering and subtraction is not equivalent to 
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the standard subtraction schemes (e.g. MS) and thus the parton distribution functions 
should in principle be modified. The modification is generally non-trivial (using e.g. 
Pythia ||] showers), however using this specific choice of showering kinematics the 
corrected (JCC) scheme is for quarks simply related to the MS parton distribution 
functions: 

zf l f I c (z,^) = zffff(z,^) (2.22) 



i 



„ ,z. , / z \ z I z 

W ? )i»(i- ? )+ ¥ (i- ¥ 



+ O (first-order quark terms) + 0{a 2 s ) . 

The above Equation |2.23| does not apply to gluons. It needs to be pointed out that 
at present it is not known how to simply extend this approach to gluons. For the 



interested reader on the work in this direction we recommend the paper [29|. 

It needs to be emphasised that for the purpose of this paper, only a very narrow imple- 
mentation of the showering algorithm was needed, namely a single heavy quark backward 
branching to gluon and subsequently only this part was actually implemented in place of 
a full parton showering program. 

We can now summarize the properties of thus obtained algorithm of combining pQCD 
and parton showering as implemented in our Monte-Carlo algorithm: 

• All heavy quarks (incoming and outgoing) are kept massive throughout the proce- 
dure, both in the perturbative calculation of matrix elements and in the showering 
procedure and overlap removal. The matrix elements used are at present 'leading 
order' (tree-level) only, however the procedure could in principle be expanded to 
diagrams containing virtual corrections. 

• All overlap removal is done at the parton level on an event-by-event basis. The 
collinear topologies are determined from the participating Feynman diagrams (cur- 
rently done manually but could in theory be automatized). 

• The kinematics of the shower and overlap removal is implemented especially for this 
approach. These choices are reflected in the subtraction term which is achieved by 
calculating the collinear limit of the kinematic topology (event) of the n th order 
perturbative calculation. 

The event generation is thus performed in the following steps: 

• The n th order process (event) is sampled, the collinear limits for the given event 
topology are estimated and the subtraction terms of order (n-1) are calculated. The 
weight is given by Equation |2.16| . 

• The (n-l) th order process (event) is sampled. If this process still contains gluon split- 
ting to heavy quarks in the initial state (in the other incoming leg) the corresponding 
subtraction terms of order (n-2) are again calculated. The event is then showered 



and the weight is given by Equation 2.13 
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• If subtraction terms were found in the previous step, the corresponding (n-2) th order 
process is calculated and showered on both legs to achieve the (n th order) event con- 
figuration of the previous two cases. The weight is estimated analogous to Equation 
p. 13 for two showering steps. 

All the above classes of events are separately unweighted, obtaining weights equal to ±1 
because in some phase space points the contributions from subtraction terms can actually 
be greater than the unsubtracted values, as will be shown in the following sections. The 



summed contributions of all processes, corresponding to Equation 2.17, are of course pos- 
itive throughout the phase space. In the actual Monte-Carlo program, it is very easy to 
(pre-)mix these classes of events internally to obtain the summed contribution correspond- 



ing to Equation 2.17 



The procedure of [14] has been improved by introducing the massive splitting kernel 
correction for the initial state g 



P, 



g-*H 



pmassless , pcorrection 



HH splits, which can be derived from [23| results: 

2z(l — z)m 2 H 



T R (1 - 2z(l - z)) + T R 



H 



(2.23) 



with pt being the relative transverse momentum of the spectator heavy quark in the split. 
The massive kernels have already been used for final state showering in the Sherpa Monte- 
Carlo generator [||] but have to our knowledge never been used in the showering of initial 
state heavy quarks. 

2.3 The gg -> Z°/~/*bb -> ffbb Implementation 

Assuming an experimentalist at LHC wants to study the Drell-Yan production with two 
resolved (heavy) b-quarks in the final state, one now has three possible choices of calculation 
(c.f. Figure ||): 




Figure 3: Schematic representation of contributions resulting in exclusive Z°HH final state: two 
fully evolved heavy (b) quarks entering 'pure' Drell-Yan at order a® in combination with double 
initial state parton shower (left), one heavy quark and one gluon entering the hard process at order 
a\ in combination with one parton shower (middle) and fully perturbative calculation involving two 
incoming gluons in a hard process of order or s (right). These three processes need to be combined 
with appropriate overlap removal as detailed in the text. 



• The order a® hard process bb — » Z°/j* — ► // with fully evolved b-quarks entering 
the hard process at \ip. The cross-section contains no mass singularities and needs 
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no subtraction terms (d&^ = da^ '). The associated b-quarks are then resolved at 
scales fjL\H and /i 2 h using parton showering: 

do AB_Z0/ 7 *66 = Y dS 9-*H{^lH) f H / f) da^^ / ^ dS g ^jj(fJ, 2 s) /h/b(Mf) • 
if=6,6 

(2.24) 

Kinematically, the leg with higher induced virtuality (scale) is treated (unresolved to 
gluon) first. 

The order a s hard process gH — > Z° /^f* —> ff H, H = b,b with one fully evolved 
b-quark entering the hard process &t fj, F and the other one participating as the propa- 
gator inside the matrix element calculation. The cross-section thus contains one mass 
singularity related to the propagator and needs a subtraction term derived from the 
collinear limit of the b-quark propagator. The other associated b-quark is then re- 
solved using parton showering and the scales \x\^ are set to be the evolution scale of 
the showered quark and the virtuality of the b-quark propagator respectively: 

da AB-*z° h*bl = Y dS 9^H{^m) Ih/a^f) da^ z0/rH f g/B {fi F ) (2.25) 

H=b,b 

+ Y fg/A(VF)da { ^ z0/rn dS g ^ R (fi 2B ) f R/B ([iF) 

H=b,b 

~ Y dSg^ H (viH) Ih/A^f) da^^ z0/r df { g' /g (fi 2lt ) fg/EiM 

H=b,b 

~ Y fa/A{^F)df^ ] /g {^ 1H )da { ^ z0/r dS g _ a (fi 2B ) /b/bM ■ 

H=b,b 

The order a 2 hard process gg — ► Z°/^/*bb — ► ffbb where both incoming b-quarks par- 
ticipate as propagators inside the matrix element calculation. The cross-section thus 
contains two mass singularities related to the propagators and needs corresponding 
subtraction terms (with permutations) derived from the collinear limit of the two 
b-quark propagators, virtualities of which then define the scales Using the 

formalism of [14|, or equivalently ||, one obtains: 



do 



2 _ _ f ..f..j\jjfi 

AB- 



zyrbb = fg/A(VF)da gg ^ z0/rb - b f g/B (fi F ) (2.26) 
- Y f9/AMdf^j g {^ 1H )da^ z0/r}1 f g/B {^ F ) 

H=b,b 

~ Y fg/A{^F)da { ^ z0h , n dff /g {ii 2R )f g / B {^ F ) 

H=b,b 

+ Y fg/A(^F)df H 1) /g (fi 1H )da < ^ l _ iZ0/r dfff /g (fi 2B ) f g/B M ■ 



H=b,b 



The derivation of the above subtraction terms is presented in Appendix [A|. 

One might be surprised at the + sign of the last 'subtraction' term in Eq. 2.26 how- 
ever re-arranging the subtraction terms to be matched with the corresponding showering 
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expressions in the spirit of Eq. 2,19 one derives the tree cross-section differentials: 
do- = ( dS 9^H(^iH) Ih/a^f) ~ df^J g (fi 1H ) fg/A^F)) da { ^^ z0/r x 



H=b,b 

< 

9- 



x [dS g _ B (fi 2a ) f B/B (fx F ) - dp R ' /g {n 2R ) f g / B (pF, 



(2.27) 



do 1 = ^2 (dSg^H^iH) Sh/a^f) - ^!s/g(PlH)fg/A(PF)jda^ 2 p /fftH f g / B (ji F ) 



H=b,b 



(1) 



.(2) 



H=b,b 
and trivially: 

do 2 = fg/A(VF)da gg ^ z0/rbi 

As noted above, there are actually four scales 
Mlif 2B1 m = ^5^) j namely the virtualities as- 
suming that the b-quark (anti-quark) are con- 
nected to the first (second) gluon from hadrons 
A and B. From the above forms it is evident 



fg/B^F) 



(2.28) 



(2.29) 



M-2 



that the first contribution from Eq. 2.27 is zero 
when either of \x\$, = fJ, F as are respective terms 
in Eq. 2.28, ensuring the smooth transitions 
of the functions in the region close to the fac- 
torization scale as expected. Subsequently, the 
behavior of the subtraction terms corresponds 
to four scenaria /Ui 2 < Hf where all the tree 
contributions are non-zero, fi± < \x F ,\i<i > fi F 
and fi2 < Mi > ^F where only the contribu- 
tions from Eqns. 2.28 and [2.29 are non-zero and 
the region > \x F where only the fully per- 
turbative contribution of Eq. 2.29| contributes. 
Graphically, the contributions are shown in Fig- 
ure |i[ 

At this point one should also review the 
omissions (limitations) of the proposed algorithm: 



a 1 + a 2 


a 2 




a 1 + a 2 



Figure 4: Schematic representation of con- 
tributions sources corresponding to four sce- 
naria fi\ t 2 < A*f where all the tree contribu- 
tions are non-zero, fi\ < > and 
H2 < RFjMi > where only the contribu- 
tions from Eq. 2.28| and 2.29| are non-zero 
and the region pb\ t 2 > where only the 



fully perturbative contribution of Eq. 2.29 
contributes. 



As presented in Equations 2.24- 2.2S) the terms dSg^nifJ-) indicate that there is no emission 



between the scales [hx^iI^f] (cf. Eq. 2.13 ), i.e. the heavy quark unresolves (backward evo- 
lution) directly to gluon. All the tree contributions will thus uniquely produce the same 
initial ( gg) and final (Z° /j* bb) state. 



In a full shower implementation as e.g. in Pythia ||, HERWIG [19] or Sherpa Q, 
one or more additional branchings in terms of gluon radiation H — ► gH could take place 
before the heavy quark would resolve back to a gluon via gluon splitting. Consequently, 
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a part of the contributions, which would be present in case of a full shower from inclu- 
sive gH — ► Z°/7* H and HH — > Z°/7* production, is missing. The probabilities of omitted 
showering contributions are however considered to be small, using again the argument of 
the order of heavy quark density fjj/i being of the effective order a s higher with respect 



to the gluon PDFs and considering the dS form of Eq. 2.13. 



3. Implementation and Results 

The process gg — > Z° /^/*bb is implemented in the AcerMC Monte-Carlo generator Q using 
the (adapted) MadGraph p4] matrix elements with full Z^/j* interference. Unweighted 



events corresponding to the three sub-processes, given by Equations 2.24, 2.25 and 2.26, are 



generated with native AcerMC single heavy quark backward branching to gluon [14|. Each 



dS term in the Eq. 2.24 , 2.25 and [2.26 thus corresponds to an actual parton showering 



step in the event generation. 

As already stated, only a very narrow implementation of the showering algorithm was 
needed for the purpose of this paper, namely a single heavy quark backward branching, and 
subsequently only this part was actually implemented in place of a full parton showering 
algorithm. The produced parton level events from AcerMC thus need to be passed to 
Monte Carlo tools as Pythia or HERWIG |l9| for further showering and hadronization. 

The showering algorithms in AcerMC (implemented in the style of Pythia veto sam- 
pling) can in principle easily be expanded to full shower with both initial and final state 
radiation and all possible branchings but for purposes of this paper we only needed that 
single showering step. Although the implementation of the full showering in AcerMC would 
guarantee consistence by using only one showering algorithm (whereas now it needs to be 
extended with further Pythia or Herwig showers with different showering algorithms ), the 
subsequent interface of the partons to the fragmentation algorithms in Pythia or Herwig 
would be rather difficult and beyond the scope of this project. 

Due to the subtraction terms a fraction of event candidates achieve negative sampling 
weights and unweighted events are produced with weight values of ±1 using the standard 
unweighting procedures. In the subsequent studies only the pure — >■ channels were 

used for benchmarking (the photon contribution was turned off) in the LHC environment 
(proton-proton collisions at y/s = 14 TeV) in combination with the CTEQ6L1 and 
hence derived JCC PDF sets 1 to manifestly see the impact of using the JCC (Eq. |2.23| ) 
modification, which turns out to be sizable in the b-quark case, contrary to its impact 
on the light quark PDF values. Let us emphasize again that Collins et al have shown 



in |TT| g0|, Ell |2^] that the PDF sets obtained e.g. in the MS subtraction scheme as the 
CTEQ sets are not formally the correct ones to be used in parton showering but that 
instead modified PDF sets corresponding to the showering algorithm (as the JCC for the 
Collins-style shower prescription implemented in AcerMC) should be used. In other words, 



1 There is a possibly valid argument that the NLO PDF sets should be used instead of LO ones since 
the derived procedure is in part NLO. The choice of these does not affect our results qualitatively and 
affects above all the absolute normalization (cross-section) prediction which is not correct anyway due to 
the absence of virtual corrections to our calculations. 
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the PDFs in Monte-Carlo event generators are determined by the showering algorithm and 
cannot be freely chosen, unlike the case for PDFs used in inclusive calculations. 

The b-quark mass is set to ruf, = 4.8 GeV and the factorization and renormalization 
scales were chosen to be equal to the Z° mass (hf = Mz) but of course other choices are 
possible. All the presented results are at parton level (i.e. b-quarks, Z /^* and its decay 
products) as output by AcerMC. 

It is instructive to first reproduce the order a), results derived in our paper |l4| in order 
to check the impact of the massive evolution kernels introduced in this paper (Eq. |2.23| ). 
The results are presented in Figure [|. 




Figure 5: The reproduced order a\ results from Q using the CTEQ6L1 and JCC PDF sets (left 
and middle); note that the subtraction term (magenta line) indeed gives a smooth transition from 
low /i region, where it matches the perturbative calculation (blue) to the /i = jip region where it 
matches the parton shower prediction (red). In the right plot the impact of introducing massive 
splitting kernels is shown explicitly for the order a\ prediction combined with the subtraction term; 
note that the low /i peak previously present in || now disappears completely. 

As one can observe in the rightmost plot of Fig. || the introduction of massive splitting 
kernels has corrected the low scale (virtuality) region; the subtraction term now indeed 
smoothly interpolates between the low scale (collinear limit) where it matches the order al 
contribution to the factorization scale [If = Mz where it coincides with the parton shower 
prediction as expected, thus allowing smooth transition between the two contributions as 
predicted. Note that at the order al two virtualities for incoming b-quarks are picked from 
the Sudakov back-evolution (to gluon) of the order a® pure Drell-Yan process but only one 
leg is actually showered and correspondingly only one subtraction term is introduced in 
the order a I process. Since in this case there is no ambiguity one can trace the actual \x 
used in the shower evolution and subtraction (labeled /i(true) in Figure ||). 

This is however not the case once we go to the full a 2 s order where, as already stated, all 
the incoming b-quarks are showered back to gluons. Here four possible scale choices exist 
Hi = \£ = —{p(gi : 2) —p(b,b)) 2 + m 2 and multiple overlap removal terms can contribute. We 
thus chose to plot the virtualities/scales related to first gluon and b-quark (/ii) and second 
gluon and anti-b quark (112) and their combinations; all other permutations in computing 
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Figure 6: The full order a 2 s results using the CTEQ6L1 (top) and JCC PDF sets (bottom), showing 
the distribution with respect to the virtuality /scale related to first gluon and b-quark (pi) (left) 
and maximal value of /ii fi2 (virtuality/scale of the second gluon and anti-b quark) (right). As 
predicted the summed contribution with overlap removal gives a smooth distribution despite the 
sharp cutoff at fx = fJ,F = Mz ■ 



the virtualities give identical predictions. The one-dimensional plots, as presented in Figure 
|(| again confirm the predictions of our approach resulting in a smooth distribution from the 
combination of different order contributions. The contributions to the total cross-sections 
are given in Table |]. 

One can observe that the contribution of the order is small both in the absolute 
normalization and its contribution to the one-dimensional projections. Its impact and 
importance is better observable plotting the two dimensional (//i,/^) differential cross- 
section plots presented in Figures ^j8[ Apart from contributing to its 'exclusive' region 
of high fj,± t 2 >> \if it contains a significant correction in the region \i\ ~ fi2 throughout 
the /ii ( 2 value range. Upon reflection this is to be expected since the perturbative-level 
calculation (and thus the full combined result) should not be biased around the [x\ ~ 
fj,2 region whereas the parton-showering approach by definition is biased in this region, 
especially at values close to fj,jr. The probability of having both flip large is the square of 
probability to obtain a single large [i value and its value is thus supposed to be low, with 
the observable dip in the distributions Figures [7|,|| which the order contribution fills 
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cr CTEQ6Ll,/i F =m z [pb] 


°JCC,/iF=m z LP 1 -*] 


a 


64.4 


44.8 


a 1 


-10.7 


-8.9 


a 2 


2.0 


2.0 




55.7 


37.9 


gg -> Zbb - 


-» n+fi-bb 22.9 


22.9 



Table 1: The process cross-sections for the leading order contribution integrated cross-section 
a , order a* contribution a 1 including the subtraction terms and order a 2 , contribution er 2 , also with 
full overlap removal, are shown for the Z° — > decay channel in the LHC environment (proton- 

proton collisions at -^/s = 14 TeV) are listed. The b-quark mass is set to mj, = 4.8 GeV and the 
factorization and renormalization scales are set to the Z° invariant mass squared. In addition, the 
perturbative (order a 2 ) gg — > Zbb — > /i + g~bb process cross-section is shown for comparison. The 
cross-sections are given for the LO CTEQ6L1 [g5| and the derived JCC PDFs. In the Monte-Carlo 
event generation procedure the next-to-lcading process weights are combined with the subtraction 
weights on the event-by-event basis as described in the text. 



up as expected. 

The kinematic quantity of interest is also the impact of subsequent corrections on the 
transverse momentum of the Z° boson and b-quarks as shown in the Figure ^. Again, the 
order a 2 contribution seems to be comparatively small but is of importance in the very 
high transverse momentum regions. As it might be, one needs the full order a 2 procedure 
to achieve the exclusive final state containing two b-quarks and formally obtain the full 
symmetry in the procedure. 



4. Conclusions and Further Perspective 



It is instructive to compare our procedure to the procedure developed in the MC@NLO 
1 26] framework. The MC@NLO Monte-Carlo generator incorporates an approach which 
results in double-counting removal both in the initial and final state — using massless 
quarks in the initial state — at full NLO, including virtual corrections, while our approach 
is currently developed to use only tree-level matrix elements but can be used iteratively 
to at least a 2 as shown in the given paper. Still, for our approach to be consistent, it 
should be compatible to the MC@NLO formalism and it will be shown this is indeed the 
case. MC@NLO is in its formalism following the paradigm of 'minimal intrusion' into the 
(showering) Monte-Carlo and is being interfaced to HERWIG Monte-Carlo generator |l9| 
with the kinematic translations adapted to its showering (evolution parameters etc . . . ). 
Fortunately, it turns out one does not need to go into details to show the consistency of 
the two approaches; it is enough to use the 'master modified subtraction equation' of the 
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Figure 7: The full order a 2 s results using the JCC PDF sets, showing the (normalized) distributions 
with respect to the virtuality/scale related to first gluon and b-quark (/xi) and virtuality/scale of the 
second gluon and anti-b quark) (/X2) for separate contributions including their subtraction terms. 
On the bottom right the sum of all contributions is shown to give a smooth distribution over the 
whole region as expected. 



MC@NLO toy model in |26|, i.e. Equation (3.20), which states: 



da_\ 
dOJ 



msub 



dx 



I M c(0,x M (x)) 



a[R(x) - BQ{x)] 



+ 7 M c(0,l) [B + aV + 



aB[Q(x) - 1] 



(4.1) 



where B represents the LO Born term, aV the finite part of the virtual corrections, aR{x)/x 
the unsubtracted real corrections and Q(x) is the Sudakov evolution kernel. The MC 'term', 



aBQ(x) 



x 



(4.2) 
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Figure 8: The full order a 2 , results using the JCC PDF sets, showing the (normalized) distributions 
with respect to the virtuality/scalc related to first gluon and b-quark (/ui) and virtuality/scale of 
the second gluon and anti-b quark) (^2) for gradual combinations of contributions of order a®' 1 ' 2 
including their subtraction terms. Note especially the disappearance of the dip in the region \i\ ~ fi2 
as detailed in the text. 



equivalent to our subtraction terms, is in the above case both subtracted and added to the 
NLO cross-section expression to obtain two separately finite contributions. As stated, in 
MC@NLO the shower evolution is assumed fixed/pre-defined by the Monte Carlo program 
(defining the x as evolution variable and Q(x) as the showering kernel) while the subtraction 
scheme at NLO, which removes the divergences and defines its own subtraction kernels is 
chosen separately (e.g. MS). For simplicity in the MC@NLO toy model the subtraction 



kernel is set simply to one and the numerator aB[Q(x) — 1] in last fraction in Eq. 4.1 



represents precisely the difference between the shower evolution kernel (scheme) and the 
NLO subtraction kernel (scheme). 



If one instead chooses another approach and adjusts the shower and NLO subtraction 



JCC_ 



do°/d(p£) <Z° -> nV) 

do^VdCpf) <z° -> n* ii) 

do 0+l+2 /d(p^)(Z°^n>-) 
do/d(p£)(g g -» Z b b -» |i + |i b b) 




■ICC 



da7d(p T w™>) (Z u - 
d0 O+l /d(p T b(m»x)) ( z" 



-» 10 



(i + n" b b) 




p» [ GeV ] 



p T b(mai) [ GeV ] 



JCC 




do"/d(p^) (Z" -» |f |l) 
dO^VdCp?) (Z° -> ^ + Hi 
do 0+l+2 /d(pZ) (Z° -» n + n") 
do/d(p^)(g g — > Z b b — > Jl + [J," b b) 



pf [ GeV ] 



Figure 9: The full order a 2 results using the JCC PDF sets, showing the differential cross-sections 
with respect to the b-quark transverse momentum (left), maximal b-quark transverse momentum 
(middle) and Z boson transverse momentum (right). The perturbative calculation of order a 2 
without subtractions is shown separately for comparison (magenta line). 



schemes to match, i.e. 



one gets instead: 



[Q(x) - 1] -» [Qshower ~ QnLo] -> 0, 



(4.3) 



da 
dO 



msub 



Imc(0,x m (x)) 



a[R(x) - BQ(x) 



+ J MC (0,1)(B + aV) 



(4.4) 



where the MC term is present only in subtraction from the real contribution but the NLO 
subtraction scheme is re-defined to match the showering scheme. If one disregards the finite 
virtual correction (which affects only the normalization) and particle masses this matches 
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exactly the formula of our approach. It also becomes clear that the parton distribution 
functions need to be re-defined since we are no longer in the MS scheme (or any other 
standard scheme) but in a shower-governed scheme as clearly explained by Collins et al 



1 21], [22 1 . 



The above comparison also shows the current limits of this approach: Our method 
assumes that there are no other divergent terms in the cross-sections apart from the mass 
divergences; furthermore, the expressions for shower scheme JCC parton distribution func- 
tions have so far been developed only for quarks. Consequently, if one wants to incorporate 
the presence of e.g. gluon radiation/splits in the hard process this approach would need 
to be extended, combined with other approaches (e.g. the addition of light 'jet-objects' 



through an algorithm like Vincia [27| or dipole showers (2£|) or possibly even take new 
directions [29|. Nevertheless, facing the stated limitations, the presented procedure is thus 
shown to be consistent with MC@NLO approach at NLO order (al corrections), while it 
does not include virtual (normalization) corrections. 



Another point of interest is to compare our approach to the one of L-CKKW [ 12 1 , [ 30 ] 
as the most advanced parton shower and matrix element matching algorithm so far. While 
L-CKKW defines a very sophisticated interpolation scheme it only picks the dominant 
collinear topology on an event-by-event basis by virtue of the employed jet clustering; in 
contrast, in our approach more than one collinearity can be accounted for and subtracted 
for each event. Furthermore, L-CKKW, while removing the double counting, was not (yet) 
shown to be formally correct in case of hadrons as the colliding particles |3l]] . On the other 
hand, the procedure described in this paper currently works only for initial state splits 
involving heavy quarks; in multi-particle final states the heavy quark lines are generally 
expected to be comparatively few. This would lead to a very complex procedure if light 
quarks/gluons were added, while it might not be necessary from the practical perspective 
since L-CKKW seems to be describing the multi- light-jet final states very well and could 
probably also be combined as a continuation of the procedure presented in this paper. 

The near-future plans are to extend this approach to the heavy quark production in the 
final state gluon splits g — ► HH and to implement the same formalism in the associated Hbb 
production which does in itself not need any further development of the formalism. There 
are certainly also plans to compare our approach with other gg — ► Z° j^*bb implementations 
and its impact to the LHC predictions is certainly envisaged but requires further work at 
the level of mock experimental LHC analysis and is thus digressing from the scope of the 
current paper. 

A. Derivation of the Hard Cross Section Expressions and Subtraction 
Terms 

The hard cross-section expressions and corresponding subtraction terms can according to 
H be actually be derived from the Factorization Theorem itself by using the Factorization 
Theorem at the parton level and doing power counting of a s . Specifically, at a given order 

(n) 

in a" the factorization Theorem implies that the perturbative pQCD cross section ^b_>X 
involving initial state partons a, b is related to its corresponding hard cross section £~kL*X 
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(containing no mass singularities and being completely infra-red safe) by the modified 
Factorization Theorem 2 : 



c,d 

where the sum of n, terms gives the correct order n, i.e. n = n\ + ri2 + n.3. The above 
equation differs from the regular Factorization Theorem expression in two respects: 

• The parton distributions are in this case relative to an on-shell parton target. 

• These parton- level distributions are expanded in powers of a s , with denoting 
the term of order a" 1 in the distribution function expansion. 

Furthermore, we limit ourselves to the case where at least one of the evolved incoming 
partons c, d is massive with mass Mh- Using the ACOT scheme [0, the above parton 
densities are calculated in the MS scheme in the region [i > Mh at the first two orders in 
a s as: 

/$«./•> -fiUfflb (A,, 

where i,j = {<?, if} and P^niO are the (usual) first order splitting kernels for (heavy) 

in) 

quark production. In addition, in this scheme, the pQCD cross-section cr ab _^ x should al- 
ready have been regularized using the MS scheme and the collinear singularities should 
have already been subtracted resulting in a finite cross section with only the mass singular- 
ities (meaning, at finite Mjj, large logarithms ln([i 2 /m H )) still present. In the calculation 
the light quark masses should be set to zero and the Mh kept at non-zero value through- 
out the calculation (contrary to the conventional approach where Mh is eventually set to 
zero). The hard cross-section o^b^x ^ s then obtained by iteratively (recursively) inverting 
the Equation [O] (c.f. 0). 

Let us illustrate how this approach was applied in this paper by explicitly considering 
the gg — > Z°/7*66 process. At order a 2 (n=2) we get from equation |A.l : 



a in^ZO/^* h T, - fala ® ^L"LxOAv*w; ® fa/a ( A ' 5 ) 



gg^Z°/r*bb J g/g gg-+Z°/i*bb J g/g 



fl=6,6 
fl=6,6 



2 In the paper || also final state evolution (i.e. fragmentation functions) are considered but in this paper 
we limit ourselves to the initial state evolution only, as stated several times in the paper. Subsequently, our 
choice of processes considered is constrained by this limitation. 
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which exhausts the possible combinations. The order n=l,2 hard cross sections o^ ' 1 ) still 
need to be determined so we move our the procedure first one order down to n=l and 
obtain (H = b,b): 



a W _ f (o) 

(0) 
H/g 

g/g 1 



+ f 



a 



a 



(l) 

Hg^Z°/7* H 

(0) 
HH- 



/, 



(0) 



g/g 
(i) 



(A.6) 



+Z /r ® fjg/gifhll) 



- = f 

gH^Z°/7*H 3 > 



<T (1) «, f (0) 

a gH^Z0/ 7 * H 59 J H/g 



Iffy 



HH^Z°/7* 



f (0) 

'ff/s' 



which still cannot be inverted so we move subsequently to n=0 and obtain: 



/ 



(o) 



(0) 



f (0) .(0) 
J H/g HH^Z°/ 7 * 



'HH^Z°/7* J H/g ^ "HH^Z°/7* 

in which case the hard and perturbative cross sections are equal, meaning: 



a 



(0) 

HH^Z°/ 7 * 



a 



(o) 

HH^Z°/7* ' 



(A.7) 



(Ai 



Inserting this result back into Equation A.6 one can now express the hard cross-section: 



Hg-»Z°/7* H 



a 



a 



(l) 

Hg^Z°/7* H 
(0) 



(A.9) 



gH^Z°/7*H 



(1) 

gH^Z°/ 7 * H 



cr 



and now these results finally into Eq. A. 5, giving: 

.(2) - „(2) 

>Z /-y*bb U gg->Z°/l*bb 

E ® P 1} 



(0) 



gg-^Z^/rfth "gg-> Z °h* 

~ E fu/g^H) ® ^Hg^Z0/ 7 *H ~ u HH^Z°/ 7 * 



(A.10) 



® / 



(i) 
H/g 



Oho) 



H=b,b 

E fH/g^H) <8 
H=b,b 

a (2) 

U gg -^Z°/-y*bb 

E « 

H=b,b 

gH^Z°/ 7 *H 

H=b,b 



(7 



(0) 
HH- 



+ Z<>/ 7 * fff/g^a) 



(A.11) 



(7 



(1) 

Hg^Z°/7* H 



fsjgi^a) 



+ E 



(7 



(0) 

HH^Z°/ 7 * 



. ® 4 1 /,(/"2ff)- 



H=b,b 
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As one can observe from Equations |A.8| , A.9| and A, 11 the hard cross section can be ex- 
pressed as the difference of the perturbative cross section and a subtraction term: 



a 



(«) 



subt ' 



(A.12) 



Let us stress again that in this approach the Mh is left at non-zero value in all the above 
terms, contrary to the conventional approach where Mh is eventually set to zero. 

It should be emphasized that in this paper only Born-level (tree-level) pQCD calcu- 
lations were used, thus implicitly assuming that there are no divergences beyond massive 
ones present for the process at hand. As a consequence, all contributions involving the 
additional participation of (soft) gluons in the final state X, which would require explicit 
NLO (loop) calculations, are not considered. 
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